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Abstract 

Boundary conditions to an existing large-eddy simulation model have been changed in 
order to simulate turbulence in the atmospheric boundary layer. Several options are now 
available, including the use of a surface energy balance. In addition, we compare 
convective boundary layer simulations with the Wangara and Minnesota field experiments 
as well as with other model results. We find excellent agreement of modelled mean profiles 
of wind and temperature with observations and good agreement for velocity variances. 
Neutral boundary layer simulation results are compared with theory and with previously 
used models. Agreement with theory is reasonable, while agreement with previous models 
is excellent. 
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1. Introduction 


Over the past two years, this group has been working in support of the numerical 
modeling arm of the NASA Wake Vortex Program. It is believed that the turbulence in the 
planetary boundary layer will have a significant effect on the evolution of wake vortices. 
Therefore, a first step is the accurate simulation of this turbulence using Large-eddy 
simulation. Eventually, a nested grid capability will enable the insertion of a wake vortex 
pair into the boundary layer. Our goal has been to add proper boundary conditions and to 
validate the TASS (Terminal Area Simulation System) model for the planetary boundary 
layer simulation. 

What follows in Section 2 is a description of the boundary condition changes to TASS 
and the associated theory. Section 3 discusses the validation results and Section 4 contains 
a discussion of neutral boundary layer runs. Section 5 is a summary and we have also 
included an appendix with instructions for using the boundary layer options with the 
model. 

2. New boundary conditions in TASS 


The horizontal velocity at the top of the computational domain may now be specified. 
This is accomplished with a three layer sponge technique much like what was done 
previously with potential temperature. This velocity may be a function of time. 

The geostrophic wind may now be specified. It is not a function of time, but may be 
a function of height. This is equivalent to specification of the horizontal pressure gradient. 

There is a choice of four possible heating (cooling) boundary conditions at the bottom 
of the computational domain: 
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1. Simple uniform heating specified as a function of time in which a rate term (in W/m^) is 
added to the equation for potential temperature at the ground. This should not be used for 
large heating rates because strong temperature gradients will not be properly accounted for 
in the surface layer similarity scheme. 

2. Specification of surface heat and moisture fluxes in kinematic units ( °K m/s and m/s). 
In this case, the rate terms are added to the equations for potential temperature and for water 
vapor at the ground. In addition, the Obukhov length is properly calculated using the value 
of the heat flux rather than using temperature gradients which may be under-resolved. 

3. Specification of the air temperature and moisture close to the ground. Surface layer 
similarity is then used to calculate the proper heat and moisture fluxes. Again, the 
Obukhov length is calculated using the surface heat flux. This method was described in 
detail in our recent Annual Report (Lin et al. 1994). 

4. Use of a surface energy budget for calculation of soil moisture, soil temperature, and 
the resulting heat and moisture fluxes to the atmosphere. We use the slab model introduced 
by Bhumralkar (1975) and Blackadar (1976). 


The details and validation of method 4 have not been shown previously and will be 
given here. In addition, we will discuss validation results for the simulation of the 
Atmospheric Boundary Layer in general. 
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2 . 1 Surface Energy Budget: Theory 


The surface energy budget is essentially equivalent to the one proposed by 
Bhumralkar (1975) and Blackadar (1976). This method was tested by Deardorff (1978) 
and shown to be both efficient and accurate. Although this addition was a central item on 
our original proposal, it was felt that the code should be as simple as possible. The first 
reason for this is that the option will be used rarely. The original plans for our research did 
not include either of the other surface boundary condition options, which are simpler and 
more accurate than a surface energy balance. To clarify, the purpose of the entire surface 
energy balance scheme is to calculate the surface heat and moisture fluxes. A large number 
of parameters describing soil characteristics goes into the scheme. Many of these are 
usually not known accurately if they are known at all. If, however, the surface heat and 
moisture flux are known with any degree of accuracy, it makes infinitely more sense to use 
those values than to try to predict them using parameters that are questionable. The same is 
true if the temperature and humidity are known at some level close to the ground. The 
second reason for keeping the scheme simple is to minimize the additional computational 
time during its use. For these reasons, a vegetation layer has not been included in the 
surface energy budget code. 

The scheme models the soil as a thin slab whose temperature changes throughout the 
day above a substrate whose temperature remains constant. The rate equation for the 
temperature of the slab can then be written as 

^ = -c,(& + L,E g - a)/(p A d,)-c 2 (7; - r„)/T, 
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where T s is the slab temperature, T m is the temperature of the substrate, c s is the heat 
capacity of the soil, T, is the diurnal period, p s is the density of the soil, d, is the depth of 

the slab, Qh is the sensible heat flux to the atmosphere, Q r is the net radiative heat flux to 
the slab, E g is the rate of evaporation of soil moisture, and L„ is the latent heat of 

evaporation, c, and c 2 are constants, set to 2 ri' 2 and 2n respectively. The last term 

represents the heat exchanged between the slab and the substrate. We may write 

= A t i ’ 

where X s is the thermal conductivity of the soil. 

The radiative heat flux may be written as 

Qr = K s + 1 l -1 t , 

where K s represents the solar radiation, iX the incoming longwave radiation, and /T the 
outgoing longwave radiation. The outward radiative flux is straightforwardly written as 

1 1 - e„&r,\ 

e g being the emmisivity of the ground and o being the Stefan-Boltzmann constant. The 

incoming longwave radiation is more complicated and must be parameterized. Like 
Deardorff, we use the approximation of Staley and Jurica (1972), 
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li={o l + a m +<J h +[\- (a, +C m +c h )][0.67][1 610q a ]° 08 }&T a 4 , 

in which q a is the mean humidity at the first model level in the atmosphere and o t , o m , and 

G h are the cloud fractions for middle, low, and high clouds, respectively. Incoming solar 
radiation is written as 


K s = (1370 — 2)^(1 -casin'? 
m 

in which a represents the albedo of the ground, Tk is the transmissivity of the atmosphere, 
and 'F is the solar zenith angle. Following Stull (1988), 

T k = (0.6 + 0.2 sin 40(1 - 0Aa h )(1 - 0.7<r ra )(1 - 0.4cx ; ) 


and 


m, 


sinT = sin0sind, -cos<j)cos8 s cos[( ^ C )-A £ ], 


where 8 s is the solar declination angle, 0 is the latitude in radians, X t is the longitude in 
radians, and t u t c is the time at 0° longitude. We may write the solar declination angle as 


* '2n(d-d r ), 

<5, = 0. cos[ — ]. 

365.25 
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where d is the day number (out of 365) and d r ( 1 73) is the day of the summer solstice. 
Q ., the sensible heat flux is estimated as 

h 

Q h = pc p u*d t . 


where 


«. = 


k47w 

ln(— )- 

Z 0 L 


and 


6 . = 


K0-e o ) 


ln(-)-^4) 

z 0 L 


In these equations, z 0 is the surface roughness length, z the height above the ground, u the 
mean eastward horizontal velocity, v the mean northward horizontal velocity, d the 

potential temperature, d B the potential temperature at z Q , and L the Obukhov length. The 

functional forms used for ^(z/L) and T w (z/L) are given in Lin et al. (1994). The problem 

which then arises is to calculate 6 0 . We use the formula of Zilitinkevich (1970), 

0 o = 0 *+^ 0 .[ O . 13(^) 0 45 ], 

k v 
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in which v is the kinematic viscosity of air, 8 S is the potential temperature of the surface 


(slab), and u, and 8 , from the previous time step are used. 

Moisture flux, Eg, is calculated in the same manner. We let 


E g = pu*q». 


where 


q. = 


fc(g-gp) 

Zo E 




and q represents the moisture fraction (g/g) in the atmosphere. 

We now come to the determination of the surface value of the moisture fraction. 
Following Deardorff (1978), a slab model is used for moisture as well as temperature. 
Thus, we let denote the soil moisture fraction in uppermost portion and w2 represent 
this value for the substrate. q s , the surface atmospheric moisture fraction is given as 


9s = a'^ a ,(^) + (l-a> a 


9s^9sa,(T s ,P s ) 


a' = min(l,w Jf !w k ), 
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where w k is the value of w at which the atmosphere is saturated and a (T P ) is the 

8 s’ sJ 

atmospheric saturation mixing ratio at the surface. The rate equation for the slab soil 
moisture is 


= (tV « ~ w 2) 


0 < w < w 

% ~ max » 


where , p w is the density of liquid water, P is the horizontally meaned precipitation rate. 

w max is the maximum va l ue of w above which runoff occurs. We let d ‘ be 10 cm and d ' 

8 1 2 

be 50 cm, typical values for most soils. We let C 2 be 0.9. C ] depends upon soil moisture 

in the following way: 


C, = 0.5 


w g /w max >0.75 


C = 14 - 22.5 [(w /w )-0. 15] 

* g max' J 


OA5<w g / Wmax <OJ5 


C. = 14 


V"W*°-15 


The rate equation for w 2 is 


dw 2 _ ~(E g - P) 
dt Pw^2 


0 < vv-) < w 

2 — rv max 
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3. Validation 


3. 1 Surface energy budget validation 

The energy budget scheme was tested as a self-standing entity (outside of TASS) 
with data from Day 1 in Lettau & Davidson (1957). This extensive field experiment was 
performed near O'Neill, Nebraska. The ground and radiative heat fluxes were measured 
directly and the latent and sensible heat fluxes could be estimated by assuming surface layer 
similarity and using the given vertical profiles of temperature and humidity. For the 
atmospheric wind speed as a function of time, we linearly interpolated between the values 
given in Lettau & Davidson, table 4.2 at a height of 6.4 m starting at 0035 local time. The 
atmospheric humidity (mixing ratio) was converted from the vapor pressure values given in 
their table 4.3.a. The linear interpolation was done in the same manner. Wind speed at the 
same height was taken from table 4.1. a and interpolated. Thus these values which are 
normally contained at the first grid level within TASS were taken from observed data for 
this validation. 

The soil moisture variables were assigned in the following way. Moisture tension 
values given in table 2.3 of Lettau & Davidson were used in conjunction with soil 
temperatures in table 2. l.a and the graph (figure 2.3.3, p. 54) to estimate the mass ratio of 
water to soil. To calculate the volume ratio, we then use 

w = ^l, 

P w 

where r m is the mass ratio, is the density of the dry soil (shown in figure 2.2. 3.2), and 

P w is the density of liquid water. We estimated values of soil moisture for 5 cm. depth 

(slab), and 40 cm. depth (substrate). The next difficulty is estimating w k . The soil at 
O’Neill, Nebraska for the referenced study was described as a "sandy loam." In Sellers 



(1965), we see that for dry clay soil, the volume ratio of solid to total is 0.417. For sand, 
the ratio is 0.585. Thus, one would estimate that for a sandy loam, the correct ratio might 
be roughly 0.5. For a sandy loam. Sellers lists the field capacity to be 0.43 times this ratio. 
The field capacity is the maximum amount of water that a soil can hold against gravity. 
Thus we have a volume of 0.22 for the field capacity. One would guess that the moisture 
fraction at which the air at the ground is saturated would be somewhat larger than the field 
capacity. Thus, we use vv=0.30. A summary of all soil parameters is shown table 1 of 
this paper. 


a 

0.25 

(m 2 /s) 

0.5 X 10' 6 

Ps c s (j/m 3 K) 

1.51 X 10 6 

£ 

"030 

W 2 

0.21 

Zs 

0.16 

Wk 1 

1X30 


Table 1 . Soil parameters used for energy budget validation case. 


Heat flux and moisture flux for the observational data were calculated from the given 
profiles by assuming surface layer similarity. In figure 1, these are plotted against the 
calculated values. Notice that, for heat flux, the observed and calculated values are nearly 
identical until about 0700, local time. After this time, the agreement is fair. Figure 2 
shows the time history of the other portions of the energy budget. All are in good 
agreement with observation with the exception of the latent heat flux. This is considerably 
larger in magnitude than the observations in the early morning and considerably smaller 
than observations in the afternoon. This discrepancy is what causes the mild disagreement 
with the sensible heat flux. Early in the morning, the larger latent heat flux results in a 
smaller sensible heat flux, and vice versa in the afternoon. This is due to the lack of a 
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Figure 1. Sensible heat flux for the validation case compared with observed values. 
Observed values were calculated from observed profiles by assuming surface layer 
similarity. 
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Figure 2. Radiative, latent, and soil heat flux for the validation case. Observed values 
were measured directly for the radiative heat fluxes, were numerically calculated in Lettau 
& Davidson for the soil heat flux from temperature profiles, and were deduced for the latent 
heat flux by assuming surface layer similarity. 


11 





vegetation parameterization in the model. Deardorff (1978) shows the same effect when 
comparing results with and without vegetation parameterization. 

Thus, as mentioned in section 1, we have chosen an energy budget scheme which is 
simple and provides acceptably accurate results. If deemed necessary in the future, we may 
add parameterization for vegetation, which would increase the accuracy. 


3.2 Wangara Validation 

Validation of TASS with data from the Wangara Experiment (Clarke et al. 1971) has 
been discussed in some detail by Lin et al. (1994). The subgrid contribution to the velocity 
variances, however, were estimated incorrectly in that paper. Thus, we will show the 
most important results here. In this case, temperature and moisture at 2m were specified. 
These values were given in the data report. Resolution for the case shown is 40X40X40, 
with a horizontal grid resolution of 125 m. and a vertical resolution of 50 m. Figure 3 
shows how well this boundary condition works for the potential temperature. Shown are 
the potential temperature profiles for 0900, 1200, and 1500 local time. The 0900 profile is 
from the observations and was used to initialize the model. Note that there is larger scatter 
in the observed data, because they represent single point measurements. The model output 
was averaged horizontally. The higher potential temperature near the surface for the 1500 
observations probably indicates that the rawinsonde balloon was released within a thermal 
plume. The mean potential temperature within the mixed layer is constant, and the 
agreement at the top of the boundary layer is excellent, meaning the balloon has probably 
moved outside of the thermal by this time. 

In figure 4, we compare the TASS results for mean winds with the observations and 
with Deardorff s (1974) results. Again, the scatter in the observed data is due to the single 
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point nature of the measurements. It is clear from the figure that TASS predicts the wind 
speeds extremely well. 



0900 observed 

— —1200 observed 

— - -1500 observed 

— - - 1200 computed 

1500 computed 


Figure 3. Observed and computed potential temperature profiles for the Wangara 
Experiment, Day 33. The 0900 profile was used to initialize the model. 


Horizontal and vertical velocity variances are shown in figure 5. The variances are 
non-dimensionalized with w*, the velocity scale for a convective boundary layer 

(w*=[g<w'&> 0 Zf<6> ] 1/3 ). Here, Z is the mixed layer depth, g is the gravitational 

acceleration, <w'&> Q [ s the heat flux at the surface, and <#> is the mean potential 

temperature within the mixed layer. Although no turbulence statistics were measured in the 
Wangara Experiment, typical dimensionless values for these variances within convective 
mixed layers are between 0.2 and 0.4, which agrees well with our results. 

In addition, we have run one simulation of the Wangara case using the energy budget 
scheme. Table 2 shows the soil parameters used. Results for potential temperature are 
shown in figure 6. Notice that the temperature specification boundary condition 
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(previously shown) is much more accurate. The energy budget results, however, are just 
as accurate as the results of Pleim & Xiu (1995), who used a similar soil model with a one 
dimensional simulation of Wangara Day 33. 

3.3 Minnesota Validation 

Because the Wangara Experiment contained no data on turbulent intensities and 
fluxes, it was necessary to look elsewhere for validation of these quantities. We chose the 
Minnesota Experiment of 1973 ( Izumi & Caughey, 1976). One of the difficulties in this 
case is that the large scale pressure gradients are not known. For example, as previously 
mentioned, the geostrophic wind profile may be used by the model. These profiles, 
however, were not measured in the experiment. To account for this forcing, we obtained 
synoptic network rawinsonde data (available every twelve hours) on the day of the 
experiment we chose to simulate. We then performed an objective analysis to extract 
geostrophic winds as a function of height within our model domain. This is extremely 
important for predicting mean horizontal winds and for comparing momentum fluxes with 
observed values. As described in Lin et al. (1994), for a steady flow, 

OZ 

/«v>-v,) = 4<kV> 
oz 

where / is the Coriolis parameter, u is the eastward velocity, v the northward 
velocity, w the vertical velocity, u g and v g denote the geostrophic components, and o 

denotes averaging. The twelve hour spaced geostrophic wind data was then interpolated in 
time to correspond to the middle of our run, remaining constant throughout the run. The 
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<U>(m/s) 



<V>(m/s) 


Figure 4. Comparison of mean horizontal velocity results from TASS with Deardorff 
(1974) and with observed data, (a) Eastward wind component and (b) northward wind 
component. Results shown are for 1200 local time, after three hours of simulation. 
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o 0.2 0.4 0.6 

<u 'u ' >/w > 2 


( b ) 



<w V '>/w 2 


Figure 5. (a) Horizontal and (b) vertical velocity variances for the TASS simulation of the 
Wangara Experiment. Values were averaged horizontally over the domain as well as over 
one hour in time, centered on the local hour indicated. Subgrid contributions are estimates 
based on the magnitude of the local deformation tensor. 
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a 

0.25 

(m 2 /s) 

0.5 X lO- 6 

PsCs (J/m 3 K) 

2.1 X lO' 6 

£ 

0.85 

w 2 

0.0245 


0.002 

W k 

0.2 


Table 2. Soil parameters used for the Wangara case 



0900 observed 

— —1200 observed 

— “ “ 1 500 observed 

— “ — 1200 computed 

1 500 computed 


<e > (K) 


Figure 6. Observed and computed potential temperature profiles for Wangara Day 33 
using the energy budget scheme in TASS. The observed 0900 profile was used to initialize 
the model. 


model used a 70X70X50 grid, with 50 meter horizontal resolution, and 36 meter vertical 
resolution with periodic horizontal boundary conditions. Flux specification was used for 
the lower boundary heating condition. 

The model was initialized with a 1 259 local time sounding from the experimental site 
on September 15, 1973. Run 5A1, with which we are comparing, contains quantities 
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averaged from 1622 to 1737 local time. Thus the model is run for over three hours before 
the comparison. Model averaging was accomplished by averaging horizontally over the 
entire domain. These averages were taken every five minutes and, in turn, averaged over 
the 75 minutes of the experiment. The experimental mixed layer height, Z ( , was 1085 
meters. The model, however, predicted a height of 1420 meters. This disparity is due 
primarily to an overestimate of the heat flux between 1259 and the observation period. 
Only the average surface heat flux during the observation period was given. 

Figure 7 shows a comparison of observed and modeled average winds. The overall 
magnitude is in good agreement, but the observed winds show a large shear within the 
mixed layer. This is most probably due to a mesoscale effect which could not be captured 
by the geostrophic wind profiles deduced from the synoptic data. This brings us to figure 
8, which shows the vertical momentum fluxes as a function of height. The flux of 

northward momentum, <v’w’>, is in excellent agreement with the observed values. The 

flux of eastward momentum, <u'w'>, is in fair agreement. The curve’s shape is similar to 
the observed profile, but the magnitudes do not agree higher up in the mixed layer. Again, 
this is a mesoscale effect and the results are quite good considering the data available for 
our synoptic forcing. 
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computed U 
— — computed V 
□ observed U 
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Figure 7. Comparison of computed mean winds with observed winds, Run 5A1 of the 
Minnesota experiment. 
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Figure 8. Same as figure 7, but for vertical momentum fluxes. 
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Figure 9 shows horizontal and vertical velocity variances. Both show good 
agreement, though there is less vertical mixing in the observed mixed layer than in the 
modeled mixed layer. This is consistent with the unusually high shear observed in figure 
7. It is normally expected that the maximum of the vertical velocity variance should occur 
at between 1/3 and 1/2 of the inversion height, as shown by the model results. The 
maximum in the observations, however, is much lower. This same effect can be seen in 
the heat flux profiles of figure 10. Here, the model shows excellent agreement low in the 
boundary layer. The observed mixed layer, however, shows an unusually small amount of 
heat flux. The modeled profile is more what one would expect given the surface heat flux 
and the inversion height. 

Thus, the Minnesota Simulation has shown TASS to be acceptable in predicting 
turbulent fluxes and variances within the convective atmospheric boundary layer. There 
seems, however, to have been some mesoscale forcing which caused lower mixing 
strengths than would be expected for a boundary layer with the given surface heat flux and 
depth. This effect could not be duplicated because the forcing was not resolved by the 
synoptic network used to obtain geostrophic wind profiles. In addition, some of the 
disagreement could be due to the lack of information of the surface heat flux as a function 
of time. Only an average was given for the experimental observation period. In the model, 
the surface heat flux time dependency was determined such that: (1) roughly the right 
amount of heat would be added from initialization to middle of the observation period to 
give the right values of potential temperature, (2) the average surface heat flux for the 
observation period would have the correct value, and (3) the time dependency would be 
consistent with diurnal variation. Other small time scale variations, however, could have 
had an effect. 
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Figure 9. (a) Horizontal and (b) vertical velocity 
Experiment, Run 5A1 . 
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Figure 10. Vertical heat flux for Run 5 A 1 of the Minnesota Experiment. 


4. Neutral boundary layer 

Although there are few high quality observational data sets of the overall structure of 
the neutral atmospheric boundary layer, we can compare with theory and empirical surface 
layer data (the surface layer is the lowest portion of the boundary layer). Given the 
postulate that stress is constant within the surface layer, it is found that the dimensionless 

wind shear, <1> M is equal to one. This wind shear is expressed as 


u* 

where u, is the friction velocity, and k is von Karmann's constant which has been 
empirically determined to be about 0.4. The friction velocity is expressed as 
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U* = Tg / P 


where x 0 is the surface stress. In addition, the dimensionless scalar gradient, <t> c should be 
one. This is expressed as 

<D _ fa <?(c) 

C C* cfe 

where -w.c, is the surface flux of the scalar c. In order to compare to other computer 
codes, we have duplicated neutral simulations performed by Andren et al. (1994), who 
compared four computer codes. 

Figure 1 1 shows profiles of O m from Andren et al. and from TASS output. We see 
that, although there is excessive shear close to the ground, TASS shows approximately the 
same results as most of those in Andren et al. The best results for O m are from Mason's 

code, which employed the backscatter method discussed in Mason and Thomson (1992). 
The thrust of this technique involves adding random velocities to the calculated velocities 
close to the ground at each time step. The rationale is that the turbulent eddies there are 
under-resolved and that this “backscatter” of energy from small scales to large scales will 
account for the lack of resolution. There are, however, a number of tunable constants in 
this process and we believe that, for simulating vortices especially, this type of method may 
lead to more problems than it solves. 

Figure 12 shows a comparison of O c between Andren et al. and TASS. Here, all 
values are quite close to one throughout the boundary layer. 
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(a) 


(b) 






Figure 11. Dimensionless wind shear profiles for a neutral boundary layer from (a) 
TASS and ( b ) Andren et al. (1994). 
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5. Summary 


In summary, we have added planetary boundary layer simulation capability to TASS 
and validated results with observed data. There are several options for the surface 
boundary conditions, one of which is a validated surface energy budget using the slab 
model. We have made comparisons with both the Wangara and Minnesota boundary layer 
experiments and have shown that results compare well with the observed data, especially 
considering the limitations in determining initial conditions. 
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Appendix : Directions for using TASS PBL boundary 


conditions 


To run TASS in the boundary layer mode, one additional input file is needed, fortran unit 
7. If this file is not present, the model will run in the original mode In the boundary layer 
mode, unit 7 must be present for all restarts, and unit 7 must contain the following 
information. 

1. Five logical variable values, format (5L4) each in the following order: UNHEAT, 
TSPEC, FLXSPEC, EBUDG, TKE. UNHEAT refers to uniform heating. If this 
variable is true, then a uniform heating is input at the first grid level. The Obukhov 
length is not explicitly calculated. TSPEC refers to temperature specification. When 
true, the heat and moisture fluxes are calculated by assuming surface layer similarity. 
The Obukhov length is explicitly calculated for stress determination. If FLXSPEC is 
true, then the heat and moisture fluxes are specified in kinematic units. Again, the 
Obukhov length is explicitly calculated. When EBUDG is true, the energy balance 
scheme is used to calculate the fluxes. Only one of the previous four variables may be 
true, otherwise an error message will appear and the run will be terminated. In 
addition, if any of the above variables are true, a random temperature perturbation is 
introduced into the first three layers of the domain to start up perturbations on a 
resolvable scale. When TKE is true, the amount of turbulent kinetic energy at each 
level in z may be specified. 

2. XITMAX, the number of data items for heating specification. 

3. Heating data. Each line is free format, with the data in the following order: 
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if UNHEAT=.T.: time in minutes, heat rate in W/m 2 , U(m/s) at top boundary, 
V(m/s) at top boundary. 

if TSPEC=.T.: time in minutes, temperature (C) at z a , humidity (g/g) at z a , 

U(m/s) at top boundary, V(m/s) at top boundary. 

if FLXSPEC=.T.: time in minutes, heat flux (K-m/s), moisture flux (m/s), U(m/s) 
at top boundary, V(m/s) at top boundary. 

if EBUDG=.T.: time in minutes, middle cloud fraction, high cloud fraction, 

U(m/s) at top boundary, V(m/s) at top boundary. 

4. If TSPEC=.T., then z a appears on the line below the heating data items. 

5. Logical variable value for GFORCE (L4). If this is true, then the geostrophic wind is 
specified and the logical variable NOSTEADY should be set to true. If GFORCE is 
false, NOSTEADY may be either true or false. 

6. If GFORCE is true, the next line must contain a logical value to be assigned to 
GWCONST (L4). If GWCONST is true, the geostrophic wind is constant with 
height and only one value of the eastward and northward components of the 
geostrophic winds need be specified. 

7. If GFORCE is true, next comes a line delimited list of geostrophic wind values, with 
the eastward component first and the northward component second on each line. 
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There must be a line for each K value starting at K=2 and ending at K=KS. If 
GWCONST is true, then only one line is necessary. 

8. If TKE =.T., the next lines contain values of turbulent kinetic energy for each K level 
starting at K=2 and ending at K=KS-2. 

9. The logical value of the variable DFLUX. If true, dust is introduced from the ground 
throughout the simulation and the environmental values throughout the atmosphere are 
0 . 

10. If DFLUX=.T., the real value of DUSTIN, the dust flux at the surface. 

1 1 . If EBUDG=.T., a line containing the values of UTC (the time at 0° longitude at the 
initialization time of the model), DAY (the day of the year, between 1 and 365), 
LNGT (the longitude), and TS (the surface temperature at initialization). 

1 2. If EBUDG=.T., a line containing the values of ALB (ground albedo) , LAMS (thermal 
conductivity of the soil in m 2 /s), TM (the substrate temperature in K), WKW (w^ in 
m 3 /m 3 ), WMAX (w max in m 3 /m 3 ), W2 (W 2 at initialization in m 3 /m 3 ), and WG (w g 
in m 3 /m 3 at initialization). 

13. If EBUDG=.T., line containing the values of CS (p*c 5 in [Jkg]/[Km 3 ]), EMISS (the 
ground emissivity), D1 PRIME ( d / in m), and D2PRIME ( ^ in m). 

We now show two examples of unit 7 files. 
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EXAMPLE 1 


In this example, KS is 13. The energy budget scheme is employed and a resolved scale 
turbulent kinetic energy is specified at each level, level 2 being the only non-zero value. 
The heating value items show no middle or high clouds. Eastward winds at the upper 
boundary vary from 0.5 m/s at the beginning of the simulation to -4.15 m/s at 999 minutes. 
Northward winds here vary from 1.10 m/s to 1.47 m/s in the same time period. There are 
1 1 heating data items. There is geostrophic wind specification which is a function of 
height. There is not dust flux at the surface and the soil parameters and other parameters 
used for the energy budget scheme appear at the bottom. 


.F. . 

F. .F. 

.T. .T. 

# unheat, tspec, flxspec, ebudg.tke 
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# xitmax: number of data items for heating 

0 . 

0.0 

0.0 0.5 1.10 


60. 

0.0 

0.0 0.07 1.96 


120. 

0.0 

0.0 -1.19 2.60 


180 

0.0 

0.0 -2.23 2.55 


240. 

0.0 

0.0 -1.67 1.96 

# time in minutes; middle cloud cover; 

300. 

0.0 

0.0 -2.85 2.10 

# high cloud cover; u at top boundary (m/s); 

360. 

0.0 

0.0 -3.13 1.89 

# v at top boundary (m/s) 

420. 

0.0 

0.0 -3.83 2.07 


480. 

0.0 

0.0 -3.74 1.79 


540. 

0.0 

0.0 -4.15 1.47 


999. 

0.0 

0.0 -4.15 1.47 


.T. 



# gforce 

.F. 



# gwconst 

-5.4275 

0.0 
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-5.2825 0.0 


# geostrophic wind values (u,v) in m/s, K=2,13 


-5.1375 0.0 

-4.9925 0.0 
-4.8475 0.0 

-4.7025 0.0 

-4.5575 0.0 

-4.4125 0.0 

-4.2675 0.0 

-4.1225 0.0 

-3.9775 0.0 

-3.8325 0.0 

0.365 # tke values 

0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
.F. 

23.34 226. -144.93 282.43 

0.25 0.5E-06 280.5 0.20 0.25 0.0245 0.002 

2.1E+06 0.85 0.10 0.50 


# dflux 

# utc, day, Ingt.ts 

# alb, lams, tm, wkw, wmax, w2, wg 

# cs, emiss, dl prime, d2prime 


EXAMPLE 2 



In this example, a neutral boundary layer is being studied. Thus the UNHEAT option 
is used with the heating rate set to 0. We specify a constant geostrophic wind of 10 m/s in 

the Eastward direction. We also introduce a dust flux of 0.001 g/gs at the lower boundary. 
Again, KS=13. 


•T. .F. .F. .F. .T. 

2 

0. 0. 10.0 0.0 
10000. 0. 10.0 0.0 
.T. 

.T. 

10.0 0.0 
0.365 
0.295 
0.245 
0.205 
0.175 
0.145 
0.120 
0.100 
0.085 
0.070 
.T. 

0.001 


# unheat, tspec, flxspec, ebudg, tke 

# xitmax: number of data items for heating 

# 0 heating rate, 10 m/s Eastward wind at top of 

# domain. 

# gforce 

# gwconst: geos, wind is constant with ht. 

# geostrophic wind 

# tke values for each vertical level 


# dflux 

# dustin 
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